Preprint typeset in JHEP style - HYPER VERSION 



CERN-PH-TH /2007-058 



o 
o 

(N 

a ■ 



Evaluating multi-loop Feynman diagrams with infrared 
and threshold singularities numerically 



(N 

(N 
> 
(N 
00 
(N 
m 
O 

o 

Oh: 



X 



Charalampos Anastasiou 

TH Unit, PH Department CERN, CH-1211 Geneva 23, Switzerland 



E-mail: 



babis@cern . ch 



Stefan Beerli and Alejandro Daleo 

Institute for Theoretical Physics, ETH, CH-8093, Zurich, Switzerland 



E-mail: 
E-mail: 



sbeerli@itp .phys . ethz . ch 



adaleo@itp . phys . ethz . ch 



Abstract: We present a method to evaluate numerically Feynman diagrams directly from 
their Feynman parameters representation. We first disentangle overlapping singularities 
using sector decomposition. Threshold singularities are treated with an appropriate contour 
deformation. We have validated our technique comparing with recent analytic results for 
the gg — > h two-loop amplitudes with heavy quarks and scalar quarks. 
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1. Introduction 



The evaluation of complicated Feynman diagrams remains one of the theoretical challenges 
to be met for the needs of the ongoing collider physics program. We develop an automated 
method to evaluate loop diagrams with infrared and threshold singularities in all kinematic 
regions. The method combines sector decomposition to simplify their infrared singular 
structures and a deformation of the integration path to treat threshold singularities. 

Sector decomposition was introduced as a simple systematic algorithm to evaluate loop 
integrals by Binoth and Heinrich J2|. The algorithm divides iteratively the integration 
region into sectors; in each sector the integration variables which could produce an over- 
lapping singularity are ordered according to their magnitude. The overlapping singularity 
takes the form of a pole in the variable which approaches the singular limit first and can 
be factored out. 

A concern for the viability of the method has been the proliferation of terms. In explicit 
calculations of cross-sections through next-to-next-to- leading order, it has been shown that 
one can write efficient sector decomposition algorithms for realistic applications in gauge 
field theories []|, [|, |5|, ||, ^, ||. However, to date, automated sector decomposition is limited 
to the calculation of infrared divergent loop diagrams in kinematic regions with trivial or 
no thresholds. This could halt progress in evaluating amplitudes for interesting scattering 
processes and fusion processes via heavy particles. 

During the last few years an inspired method is being developed for the evaluation of 
one-loop amplitudes by Nagy and Soper ||, [10|] . In their method, the infrared and ultra- 
violet divergences are matched algorithmically by simple counterterms for each diagram; 
these add up to functions which integrate to the known universal poles in e of one-loop 
amplitudes. After the one-loop integrals are rendered finite in four dimensions, they per- 
form a numerical integration over Feynman parameters and the loop momentum. They 
have proposed a systematic way to find an integration contour in the space of Feynman 
parameters which is suitable for a numerical integration. 

We have merged the algorithm for sector decomposition with the contour deformation 
of Feynman parameters proposed by Nagy and Soper . Since sector decomposition offers 
a general solution to rendering Feynman diagrams with divergences in d — > 4 dimensions 
finite, in principle, we can now compute generic multi-loop integrals numerically in all 
kinematic regions. 

We have written three independent computer implementations of our method and per- 
formed extensive checks evaluating a variety of one and two-loop scalar and tensor integrals 
which we could verify with other methods. To prove the efficiency of our method, we have 
recomputed all diagrams in the two-loop amplitudes for gg — > h production via heavy 
quarks and squarks, recently computed analytically. The new method yields numerical 
results which are in excellent agreement with the analytic evaluation |TT|. 



Lazopoulos, Melnikov and Petriello have recently presented [12| the evaluation of the 
NLO QCD corrections to pp — > ZZZ. The method in their publication is the same as the 
one we are presenting. Here, we are applying it to the evaluation of a two-loop amplitude; 
together with the work of [12] this emphasizes further the versatility of the method. 
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Binoth et al. have used a contour deformation to evaluate integrals which are free of 



infrared divergences [13|. They also noted the possibility of merging sector decomposition 
and contour deformation as an alternative to their reduction method. To our understand- 
ing, the viability of this idea for practical applications was not investigated in [|D|]. 

A competitive numerical method for the evaluation of loop amplitudes is via Mellin- 
Barnes representations [14, Loop integrals with a small number of kinematic scales tend 



to have Mellin-Barnes representations with a lower dimensionality than the corresponding 
Feynman parameter representations. In such cases the Mellin-Barnes method should be 
advantageous; however, the method of this paper could perform better by increasing the 



number of kinematic scales. In addition, as noted in [15], Mellin-Barnes integrals cannot 
be computed stochastically in phase-space regions with mass dependent thresholds. For 
such applications the Mellin-Barnes method can be used only for checking purposes in the 
Euclidean region; sector decomposition with contour deformation could be then the only 
viable numerical method to obtain a physical result. 

We should also note that a very significant progress in developing numerically inte- 
grable representations for two-loop three-point functions with threshold and infrared sin- 
gularities has been made in |l6| . However, this approach is yet not fully automated or 
general and requires an indivindual study for each topology. 

We now present the method and our numerical results and comparisons. 

2. Method 

We first introduce independent Feynman parameters for a multi-loop Feynman diagram. 
In general, this yields a sum of terms of the form 



5 ^°Jo ' ~~"[g@,M?,8 M )-i8]' 



C(e)hm/ d Xl ---dx n n J \' > lc (2.1) 



where a is an integer, ni is the number of loops and Mi,Ski are masses and kinematic 
invariants. The function Q is a polynomial of the independent Feynman parameters x. 

The integrand can be singular at the edges of the integration region. As a first step, we 
disentangle overlapping singularities using sector decomposition [17, 18, [jj]. The outcome 
is a sum of integrals of the type 

' * * dXfiX-^ ' ' ' Xy\ s ; ^0 

'^Jo [g s {x,Mf,s kl )-i8\ 



I s = C(e) lrm / T^TTTo ; .^a+n L e ( 2 - 2 ) 



The singularities at the edges of the integration region are now factorized. The function 
Q s is finite at Xi — ► 0, 1 and the singularities from the factors x~ a+ @ e when a > can be 
extracted independently for each integration variable. However, Q s may produce singular- 
ities if it vanishes inside the integration region. It is important to postpone the extraction 
of the e poles until we treat these threshold singularities first. 

Following the method of Nagy and Soper (]Tc| ], we construct a contour of integration 
where the imaginary part of Q s is negative, that is, enforcing the — % 5 prescription already 
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present in the original integral. The new contour is defined by deforming the integration 
path of every Feynman parameter; it is parameterized by 

dQ s 

Zi = Xi - i\xi(l - x i)-g^ ■ ( 2 -3) 
Provided that no poles were crossed in going from [0, 1] to the contour C, we have 
-i (n? =1 dxjxp + ^) F,(2, e) f (U] =1 dzjzp' *') F s (z, e) 

(2.4) 

The choice of Eq. (|2.3|) for the contour deformation guarantees that for small values of A, 
the function Q s acquires a negative imaginary part of order 0(A) 

g s (z) = g s (x) xft - Xl ) 2 + o(\ 2 ), (2.5) 

where the 0(A 2 ) terms are purely real. One could add higher order A™ terms in Eq. (|2.5| ) 
to cancel the imaginary parts of Q s at C(A 3 ) and higher orders. In practice, it is sufficient 
to perform a linear deformation as in Eq. (^^), and choose a small enough value for A such 
that these contributions are suppressed. 

Changing variables using the parameterization in Eq. ( |2.3| ) , each sector integral is now 
written as 

„! n 



Is = C(e) / Hdx jZ p + ^ J(x ^ z)C(z(x),e) 

J ■ -i 



C(e) / J] te&T - J{* ~+ e) (2.6) 



o 



where J{x — > z) is the Jacobian of the transformation of Eq. (2.3). The function £ is finite 
in the boundary of the integration region and thus can be expanded as a Taylor series 
around e = 0. From Eq. (|2,3p we can see that the ratios Zi/xi = 1 + 0(xi) are also smooth 
in the singular limits. 

Now, we are free to extract the e poles in each integration variable by applying 



1 dxx- n+ * f{x) = ! l dxx J{x) gLo^JZ + y / (fc) (0) (2 ?) 

y ' Jo x n £^ o k\{k + l-n + e) V ' 



on Eq. ( |2.6[ ). After this expansion, we are left with integrals which can be safely expanded 
in power series in e. We compute the coefficients of the e series numerically. Notice that in 
presence of higher order singularities, n > 1, the expansion in Eq. Q2.7| ) involves derivatives 
of both the Jacobian and the function C that contains all non-singular factors, including 
factors coming from the tensor structure of the integral. 

There are many options for the numerical evaluation of the resulting integrals. For 
example, we can combine the contributions from all sectors into a single integrand; alterna- 
tively, we can integrate each sector separately and sum up the results. We have found that 
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the second choice is usually better, since the adaptation of numerical integration algorithms 
is more effective when dealing with simpler integrands. 

We have implemented the method described above into three different programs. Sec- 
tor decomposition and contour deformation are performed with MAPLE and MATHE- 
MATICA routines. The same programs control the creation of numerical routines for the 
evaluation of the integrals; these are written in FORTRAN or C++. In all implementa- 
tions we use the integration routines in the Cuba [19] library, relying mostly on the Cuhre 
and Divonne algorithms. 

An important issue when performing the numerical integration is the value of the 
parameter A in Eq. ( |2.3| ), which controls the magnitude of the contour deformation. A 
very small value could result to instabilities due to rounding errors. A very big value could 
result to a deformation with the wrong sign in Eq. ( |2.3| ) . This last case is easy to detect at 
runtime and we have implemented diagnostic routines to abort the numerical evaluation if 
Q s is found to have an imaginary part with the wrong sign. As we will show in the next 
section in an specific example, there is usually a very comfortable interval for A where the 
result of the integration is insensitive to its value. 

3. Results 

We have applied our method to compute the two-loop amplitudes for gg — > h mediated 
by a heavy quark or a scalar quark purely numerically. These amplitudes have been 



computed earlier either by using a mixture of analytic and numerical integrations [2C, 21] 



or analytically in |]Tl[ and in [22] 1 . We have evaluated all Feynman diagrams in these 





(a) (b) 

Figure 1: (a) Feynman diagram contributing to gg — > h with a heavy quark loop, (b) Master 
integral arising in the calculation of gg — > h in the MSSM. 

amplitudes with a very good numerical precision. As an example, in Fig. || we show our 
results for the diagram in Fig. |l|a. On the left panel we plot the real part of the finite 
piece of the diagram as a function of r = s/(4m 2 ), normalized to m 2 = 1. The inset 
plot gives a more detailed view of the threshold region, where the numerical integration is 



most difficult, superimposed with the analytic results from [11|. The plot on the right panel 
shows the percent difference between the numerical result and the analytic one, normalized 
to the analytic value. In black lines we included the bands corresponding to the integration 
error, obtained by adding in quadrature the errors quoted by the integration routine for 



lr rhe integral representation of the result in Ji^] was expressed in an analytic form in 
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Figure 2: Results for the real part of the finite piece (Re(co)) of the Feynman diagram in Fig. [l]a. 
The left panel shows the results of the numerical integration as black dots with error bars. The inset 
plot zooms in on the threshold region, the red line corresponds to the evaluation of the analytic 
result of fll|| . The right panel shows the difference in percent of the numerical evaluation and 
the analytic one, normalized to the latter. The gray bands correspond to the integration error. At 
threshold this error is 3%. 



each sector. We obtain similar results for the single pole and for the imaginary parts of 
both the single pole and the finite piece of this diagram. 

As stressed above, the method relies on the proper choice of the value for the parameter 
A. Values that are too large produce imaginary parts with the wrong sign for the function 
Q s . Very small values generate a contour that is too close to the real line -and thus to 
the zeros of Q s - and produce numerical instabilities. In practical implementations, we 
found that there is usually a good range for A. As an example, in Fig. |3| we plot the 
results for the scalar integral corresponding to the diagram in Fig. [l|a as a function of A 
for two different values of r. The results in these plots show that away from the threshold 
region, the integration is rather insensitive to the value of A chosen, as long as it induces 
a deformation providing the right sign for the imaginary part in Eq. (2.5). On the other 



hand, close to the threshold region, the magnitude of the deformation has to be larger in 
order to get a reliable estimate of the integral. 

As a novel result, we applied our method to the calculation of the scalar integral 
corresponding to the Feynman diagram of Fig. |l]b. This integral is one of the master 
integrals appearing in the SUSY QCD corrections to gg — ► h and it involves a massive 
quark, a massive scalar quark and a massive gluino. Our results are displayed in Fig. || as 
a function of r = s/(4m^) for fixed values of m| = 400/175mg and m~ = 600/175 with 
m q = 1. The results are again very stable over a wide range of A. Due to the absence of 
massless propagators, the numerical evaluation of this integral turns out to be substantially 
faster than the scalar integral corresponding to the diagram of Fig. @a. 
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Figure 3: Results for the real part of the finite piece of the scalar integral corresponding to the 
Feynman diagram in Fig. |l|a as a function of the parameter A for two fixed values of the kinematical 
ratio r. The results have been normalized to the analytic result for this master integral obtained 
in 0. 
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Figure 4: Results for the scalar integral corresponding to the Feynman diagram in Fig. |l|b as a 
function of t = s/(Am q ) for fixed values of m| = 400/175 m 2 q and m? = 600/175 rr? q . The inset 
plots zoom in the threshold region. The estimated relative accuracy of the points is better that 1 
per millc. 



4. Conclusions 

We present a new method for the numerical evaluation of multi-loop Feynman diagrams 
containing both infrared and threshold singularities. The method uses sector decomposi- 
tion to extract the infrared singularities followed by contour deformation in the Feynman 
parameters to deal with the thresholds present in the diagram. The algorithmic nature of 
the approach naturally leads to a high degree of automatization in all the stages of the 
calculation. 

We tested the method recalculating the two loop corrections to gg — > h mediated by a 
massive quark and a massive scalar quark. We find that the method is very efficient and 
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reliable, reproducing the analytic results with great accuracy. 

Currently we are applying the technique to the calculation of the two loop SUSY QCD 
corrections to gg — * h. As an example, we have presented here results for one of the 
most complicated with analytical methods, yet uncalculated, master integrals appearing in 
this amplitude. We find excellent numerical behavior, showing that the framework has a 
great potential for computing automatically general multi-loop processes involving internal 
thresholds. 
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